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The vertical trajectory optimization for the en route descent phase is studied in the pres¬ 
ence of both along track and cross winds, which are both modeled as functions of altitude. 
The flight range covers some portion of a cruise segment and ends at a meter flx. The 
descent trajectory is assumed as a flight idle descent. The problem is formulated as an op¬ 
timal control problem with both mixed state path constraints and pure state constraints. 
For minimizing environmental impacts, we optimize descent trajectory with respect to two 
cost functionals: fuel and emissions cost. We analyze both singular arc and boundary arc 
using the necessary conditions of the optimality. From the analysis result, we propose the 
optimal trajectory generation method, in which the optimal trajectory is generated by the 
backward and forward integration. The trajectory from the proposed algorithm is com¬ 
pared to the numerical optimal solution of the original optimal control problem. The result 
shows that proposed algorithm generates the same solution as the numerical optimal so¬ 
lution. The wind speed, wind shear, and cross wind effects on the optimal trajectory are 
analyzed with two aircraft type, Boeing 737-500 and Boeing 767-400. 
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Nomenclature 

Drag force 
Altitude 
Cruise cost 
Idle descent cost 
Lift force 
Mach number 
Aircraft mass 
Thrust force 
Calibrated airspeed 
Ground speed 
True airspeed 

Along track component of the wind velocity 
Cross track component of the wind velocity 
Along track distance 
Vector expressed in relative wind frame 
Aerodynamic flight path angle 
Heading angle 

Angle between airspeed and ground speed 
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I. Introduction 


T he Next Generation Air Transportation System (NextGen) is being developed and implemented to 
enhance the efficiency and safety of air transportation and reduce the environmental impact due to noise 
and aircraft engine gaseous emissions. The Continuous Descent Arrival (CDA) or Optimized Profile Descent 
(OPD) is a promising approach to reduce the environmental impact of aircraft operations[l, 2, 3]. During 
CDA\OPD, aircraft descend from the cruise altitude to the Final Approach Fix (FAF) at or near idle thrust 
without level segments at low altitude, thereby minimizing the need for high thrust levels to remain at a 
constant altitude, and reducing the environmental impact[4]. 

In order to maximize the environmental benefits of CDA, Park and Clarke [5, 6] formulated the vertical 
trajectory optimization problem as a multi-phase optimal control problem and showed that the possible 
trajectory variation with respect to the various performance indices occurs in en route descent phase, which 
starts during the latter stages of the cruise segment and ends at the meter fix. Thus, the focus of this paper is 
the vertical trajectory generation problem in the en route descent phase. 

Optimal control techniques have been widely used to solve trajectory generation problems. In 1980s, 
most approaches used energy state approximation in which the energy state was defined by the combination 
of height and speed to reduce the problem size[7, 8, 9, 10, 1 1]. Erzberger and Fee[7] proposed an optimal 
trajectory generation method with specified range. They used Direct Operational Cost (DOC) combining 
time cost and fuel cost as a performance index. Sorenson and Waters [8] addressed the fuel optimal trajectory 
problem with fixed arrival time. Available time delay analysis was performed by setting negative time cost. 
Burrows[9] converted the fixed arrival time problem to a DOC optimal trajectory problem with free final 
fime. Using fhis approach, the fuel optimal control problem is equivalent to finding a time cost for which 
a corresponding free final fime DOC opfimal trajectory arrives at the assigned time. Chakravarty[10] used 
the singular perturbation method as well as an energy state approximation and investigated the sensitivity of 
the fuel optimal trajectories to wind. Shultz[l 1] studied the three dimensional minimum time problem with 
fixed initial and final points. 

Energy state approximation is beneficial because it enables quicker computation and thus lower solution 
times. However, since the model is simplified, and some consfraints in the real procedure are not considered, 
trajectories generated by these approaches are typically not accurate enough to apply directly in real air 
traffic situations. 

From early 2000s, several studies for trajectory optimization that minimize environmental impact have 
been conducted. Visser and Wijnen [12] addressed noise abatement trajectory optimization problem. They 
formulated this problem as a multi-phase optimal control problem and solved by using direct numerical 
method. Wu and Zhao [13] assumed descent trajectory with several segments and converted optimal control 
problem to parameter optimization problems with fuel and emission costs. Few papers used CDA trajectory 
structure assuming idle thrust descent [14, 5, 15]. Franco et al.[14] solved maximum descent range problem 
by singular arc analysis in the presence of the along track wind. However, they did not consider the cross 
wind and the flight envelope which is modeled as a pure state inequality constraint. Zhao and Tsiotras [16] 
proposed optimal speed profile generation method for minimizing energy with given path using singular arc 
analysis. However, they did not consider the wind effect. 

This paper considers the optimal trajectory that minimizes environmental impact such as fuel burn and 
emissions. The flight range considered in this paper is also fixed from initial point in the cruise segment 
to the meter fix determined by Standard Terminal Arrival Route (STAR) procedure. In the previous work 
from same authors, optimal control framework for vertical trajectory optimization with various performance 
indices was developed[5], and fuel optimal trajectory generation method using Flight Management System 
(FMS) Vertical Navigation (VNAV) function was developed and solved in the hybrid system framework[17]. 
In this paper, we analyze the optimal solutions with respect to the fuel and emission performance indices by 
using the necessary conditions of the optimal control problem with state inequality path constraints. 

The remaining part of the paper is structured as follows: The CDA trajectory structure and optimal 
control problem are presented in section II. The analysis of the optimal solution structure including the 
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singular arc and boundary arc are presented in section III. The fast trajectory generation method is presented 
in section IV. Numerical examples with various wind conditions are presented in section V. The conclusions 
of this study are presented in section VI. 

II. Optimal Control Problem Formulation 
A. CDA trajectory structure 

The CDA trajectory structure is shown in Fig. 1 . As shown in the figure, the CDA procedure begins towards 
the end of the cruise segment and ends at either the FAF or the Initial Approach Fix (lAF). The traffic flow is 
mefered af fhe enfry poinfs of fhe terminal area, where speed and altitude constraints are typically prescribed. 
Therefore, the CDA trajectory can be divided into two parts: an en route descent phase and a terminal area 
descent and approach. Only the en route descent is considered in this paper. It consists of two segments: 
cruise and descent. Typically, the lateral path is defined via a STAR fhaf provides fhe waypoinfs for fhe 
Laferal Navigafion (LNAV) function in fhe FMS [18]. Thus, only fhe verfical pafh is free fo be optimized. 



Figure 1. CDA trajectory structure 


B. Flight Dynamic Model 

The relation between airspeed, cross wind, and along track wind is defined in Fig. 2. The angles ijj and 'ipw 
denote the heading angle and the angle between airspeed and ground speed, respectively. 

As previously stated, only vertical motion need be considered for trajectory optimization because the 
lateral path of an arrival procedure is typically prescribed. In addition, if the flight path angle is bounded 
by small values between — 6 ° and 0 °, sin 7 and cosy can be approximated as 7 and 1 , respectively, and the 
equations of motion in the vertical plane in [5] can be simplified as follows: 

VT = ^iT-D)-gj- Vt7 [c{Vt, h)^ + s{Vt, (1) 

X, = c{Vt, h)VT + Wh (2) 

h = Vt 7 (3) 

where, Vt is the true airspeed of aircraft, Xs is the along track distance from runway threshold; h is the 
altitude; 7 is the aerodynamic flight path angle; Wh and Wc are the along track wind speed and cross wind 
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Figure 2. Relation between airspeed, ground speed, and wind vector 


speed, respectively; dWh/dh and dWc/dh are the wind shear terms, c{Vt, h) = cosif^w, and s{Vt, h) = 
sin Both the along track wind and the cross wind are assumed to be functions of altitude. The aircraft 
mass m is assumed to be a constant during the entire en route descent phase. The lift force is assumed to be 
equal to the gravity force. The simplified lift force L with given assumption and drag force D are as follows: 


= mg 

D{h,VT) = ^p{h)VSSCD{VT,h,CL) 


( 4 ) 

( 5 ) 


where p is the air density which is the function of h, and Cl and Cd are the lift and drag coefficients, 
respectively. The International Standard Atmosphere (ISA) model is used. From Eq. (4), Cl = Cl{Vt, h). 
Hence, the drag in Eq. (5) is a function of Vt and h also. 

During the descent segment, thrust T is set to idle as shown in Eig. 1. Generally, idle thrust is mod¬ 
eled as a function of Vt and h using data from aircraft performance models such as the performance en¬ 
gineering manual provided to airlines by the manufacturers and the Base of Aircraft Data (BADA) from 
Eurocontrol[19]. Therefore, 7 is the only control input during the idle descent segment. 


C. Environmental Cost Indices 

In this study, fuel bum cost and emissions cost are used as environmental cost indices. The cost functional 
is defined as a summation of the cost of the cruise segment and that of the idle descent segment. Erom the 
wind speed assumption, the ground speed in the cruise phase is a constant. Therefore, the cost of the cruise 
phase can be expressed as a Mayer cost term. Eor this reason, the optimal control problem for en route 
descent part of CDA is formulated as a single phase optimal control problem with a cost functional of the 
form 

J — Kcr(xsito) dmax) T Kdg^giVT 1 h) dt (6) 

to 

where K^r is the cruise cost coefficient which is the cost per distance, and K^es is the idle descent cost 
coefficient which is the cost per time; dmax is the along track distance from runway to the initial point; 
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Xs{to) is the Top Of Descent (TOD) point in Fig. 1. The fuel cost and emissions cost can be expressed as 
follows: 


■ Fuel burn cost 


■ Emissions cost 


Kdes{VT,h) = fidleiVTM 


j. EIxfcr 
Kcr = —^ - 

^ cr 

KdesiVr, h) = EIx{Vt, h)Adie{VT, h) 


( 7 ) 


( 8 ) 


where fcr and fidie denote fuel flow rate in the cruise segment and the idle descent segment, respectively; 
EIx denotes emission index, and subscript X denotes emission gas types: CO, HC, and NOx- 

The emissions index table can be obtained from the International Civil Aviation Organization (ICAO) 
engine exhaust emissions databank[20], which is subsequently interpolated using Boeing method2 (BM2), 
an El correction method considering airspeed and atmosphere conditions such as temperature, pressure, 
and humidity[21]. Since the idle fuel flow rate is a function of Mach number and h, El obtained by BM2 
is also a function of Vt and h. Therefore, Kdes in Eq. (6) can be expressed as Kdes{VT, h) in both fuel and 
emission cost cases. 


D. Constraints for CDA 

We consider several path constraints for flight envelope protection and passenger comfort. Bounds are 
placed on the Calibrated airspeed (CAS) and the Mach number for flight envelope protection. CAS and 
Mach are functions of Vt and h as below: 


VcAsiVT,h) 


7R{Qq)isa 


(^1 + P 5 {h) 



7Re{h) 



0.5 


( 9 ) 


M{Vt, h) = Er/Vl-4i?0(/i) (10) 

where Ps is a pressure ratio between current altitude and mean sea level (MSL); i? is a gas constant; (0o)/5A 
and 0 denote the standard atmospheric temperature at MSL and the temperature at the current altitude, 
respectively. Hence, Mach/CAS bounds are pure state inequality constraints. 

For passenger comfort, we add a Rate of Descent (ROD) bound and a flight path angle bound. Boundary 
conditions are given in Eq. (15) and (16). CAS and altitude constraints are usually given at meter Axes 
per the published STAR chart. CAS constraint with an altitude constraint can be converted to the fixed Vt 
condition at the meter fix. 


■ Elighf Envelope Profecfion 


ymin,CAS < VcAS < Vmax,CAS 

Mmin <M< Umax 


( 11 ) 

( 12 ) 


■ Passenger Comforf 


RODrain < < RODmax 

at 


'Jrain ^ T — 'Imax 


(13) 

(14) 
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■ Boundary conditions 


(15) 

(16) 


Vrito) = Vto, h{to) = ho, Xs{to) = free, fo = 0 
Vritf) = Vrf, h{tf) = hf, Xsitf) = Sf, tj = free 


III. Analysis of the Optimal Trajectory 

A. Equivalent Optimal Control Problem 

Consider the following optimal control problem: 

Problem 1. 

min •^ = y -^cricVr + Wh) + KdesiVr, h) dt 
to 

s.t. (17) 

dynamic constraints: Eq. (1) (3) 

path constraints : Eq. (11) rsj (14) 
boundary conditions : Eq. (15) and (16) 

The following lemma shows that Problem 1 is equivalent to the original problem with Eq. ( 6 ). 

Lemma 1. Let x*(fo) be the optimal TOD of Problem 1. If x*{to) < dmax, Problem 1 is equivalent to the 

original optimal problem with performance index in Eq. (6). 

Proof Erom Eq. (2), Xgfo) can be expressed as 


Xs{to) 



cVt + lE/j dt. 


to 


Substituting Eq. (18) into Eq. ( 6 ), then 


(18) 


*/ 

J = KcriXsitf) - dmax) + j -KcricVr + Wh) + KdesiVr, h) dt (19) 

*0 

Since Kcrixsitf) — dmax) is a constant from the terminal condition in Eq. (16), Problem 1 is equivalent to 
the original problem with Eq. (6). □ 

Note that in the last part of this paper, we use Problem 1 to derive the necessary conditions of the original 
problem with cost in (6) from lemma 1 . 

B. Necessary Conditions 

The Hamiltonian for Problem 1 is defined as follows: 

H = -KericVr + Wh) + KdesiVr, h) + \v{^^ - 57 - VriVTh,^) + KicVr + Wh) + XhVrl (20) 

m 

where Ay, A^,, and \h are the costate variables corresponding to the state variables Vr, Xg and h, respectively. 
In Eq. (20), VVh^x ~ 3" paper, we define fhe Eagrangian using fhe direcf adjoining 

approach in [22] . The Eagrangian L including pafh consfrainfs is given by 

L = H + p^CiVr,h,-i) + q^SiVr,h) ( 21 ) 
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where the inequality eonstraints C{Vt, h, 7 ) : —;■ M”? and S{Vt, /i) : ^ represent the mixed path 
constraints and the pure state inequality in Problem 1, respectively; ^ G and 77 G are the Lagrange 
multiplier for mixed and pure state inequality constraints, respectively. From the constraints in Eqs. (11) ~ 
(14), C{Vt, h, 7 ) and S{Vt, h) have four components as follows: 


C = 


'-h{VT,j)- RODmax' 
RODmin + h(VT, 7 ) 

<0 S' — 

VcAsiVT, h) - Vmax,CAS 
Vmin,CAS “ VcAsiVT, h) 

T 'Imax 


M{Vt, h) — Mmax 

'Jmin T 


Mmin — M{Vt, h) 


< 0 


( 22 ) 


The following equations are the necessary conditions for the optimality of the Problem 1[23]. 
From the Euler-Lagrange equations. 


Ay = - 


dL 


dVr 


dc 


dK, 


= {K„ - y) I c + I - ^ + M 


1 dD 


dWh- 


dC 




dVr 


d\ 


t dL ^ 
Ax — — ^— — 0 
UXg 


(23) 

(24) 


dL 




dc. 


dWh\ dKdes 


1 dD 


dh dh 


dWh,x^ „TdC jrdS 




(25) 


where D = D — T, which means the net drag force. 

Erom the Karush-Kuhn-Tucker (KKT) conditions, 

fi > 0, i/'C{VT,h,'y) = 0, 

r?>0, r]^S{VT,h) = 0. 

The optimal control is determined by 

7° = argmin^gf^(y^ 

L'y = H'y = 0 


(26) 


(27) 

(28) 


where Q{Vt, /i) = { 7 | C{Vt, (i, 7 ) < 0 } is the admissible control set, and it depends on the state variable 
Vt and h. Since Xs{to) and final time fy are free, the following transversality conditions hold: 

Ax(fo) = 0 (29) 

H{tf) = 0 (30) 


Since state inequalities are independent to time and Xg, for any time r on the state boundary arc where 
S{Vt, h) = 0 , the following jump conditions hold[23]: 

Ax(t") = A 3 ;(t+) (31) 

H{t-) = H{t+) (32) 


Therefore A^; and H are continuous function with respect to time. Combining with Eq. (24), this condition 
leads to 

Aa,(t) = 0 for t€[to,tf]. (33) 

Eurthermore, since the Hamiltonian in Eq. (21) is not an explicit function of time and the cost index in 
Eq. (17) is in a Eagrange form, the Hamiltonian is a constant along the optimal trajectory. Hence, from 
Eq. (30), 

H{t) = 0 for t£[to,tf] (34) 
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However, Ay and Xh can be discontinuous for any time r in the boundary arc interval and for any contact 
time r by the following jump conditions [23]: 


(99 

Ay(r ) = Ay(r+) gy^ 

(35) 

A/i(r") = A,j(t+) - y(r)^|^ 

(36) 

y(r) > 0 , y(r)5(r) = 0 . 

(37) 


To solve this problem with the necessary conditions in Eq. (23)~ (34), both the interior arc on which 
S{Vt, h) < 0 and the boundary arc on which at least one of the components of S{Vt, h) is zero must be 
considered [24]. 


C. Singular Arc 

First, assume that all the state constraints S{Vt, h) are not active. Then rj is set to zero. By the minimum 
principle in Eq. (27), the optimal control input 7 is determined as 


Tmax 

A 

0 


7 

if 77.^ = 0 

(38) 

7min 

if 77.^ > 0 



where 7 is the singular control, and 7max and 7min are determined by 0,{Vt, h). The switching function 
in Eq. (38) is given by 

= -Xvig + VtWh,^) + XhVr. (39) 

For the analysis of the singular arc, assume that H,y is zero during a finite time interval. Then the singular 
control and the singular arc are obtained by the time derivatives of From = 0, we can obtain the 
following relation: 

{g + VTWh,^) = ^VT (40) 

Ay 

By differentiating Eq. (39), and substituting for the time derivatives of the costate and state variables 
and Eq. (40), is expressed as 


H., = Kc 


- c + 


W/t) 


m 


D Wh,x + 


dW 


h,X 


dVr 


dh 


dKdes dKdes-rr 

H ^— Vt - TT, — Vt 


dVr Ay 


dh 



(41) 


- — -TwryT + D] = 0 . 


From Eq. (34) and the singular control assumption H-y = 0, the remaining term in Eq. (20) is zero on 
the singular arc also. Therefore, 


Ho = -KcricVT + Wh) + Kdes 


, D 

Ay— = 0 
m 


By combining Eq. (42) and Eq. (40), 


Ay 

m 

bh. 

m 


i{-KcricVT + Wh) + Kdes) 
m Vt 


(42) 


(43) 
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By substituting Eq. (43) into Eq. (41), the final formulation of the singular arc rs(Vr, h) is obtained as 


FsiVT, h) = [K,r{cVT + Wh) 


KdesiVrM 



db 

dVj' 


(g + VT^h,x) + D 


d^h,x 

dVr 


- D 

- b 
= 0 


(cs ^Vrj (g + if + FrW,,,) 

Kdes dKdes / -TIT \ dKdes^T 


(44) 

Since the singular arc rs(V 7 ’, h) in Eq. (44) is not an explicit function of control 7 , second time derivative 
is needed to obtain the analytic formulation of optimal singular control 7 . Eurthermore, the following 
Generalized Eegendre-Clebsch (GEC) condition should hold if the singular arc is a part of the optimal 
trajectory. 


A 

07 



A 

d'j 

A 

d'j 


dt 

dH^ 

Wr 


7 


V) \ r)M 

---(p + EtW7,)7 +^(Et7) 


+ VTWh,x) + ^Vt <0 


dVr 


dh 


(45) 


Here, the partial derivative terms of can be calculated numerically when Eq. (44) is solved. Hence, this 
condition can be checked on the singular arc. All numerical examples in this paper satisfy this condition. 
Eurthermore, if the singular control holds 7 min < 7 < 7 max> the singular arc control satisfies the following 
junction theorem also. 


Theorem 1 (Junction Thm in [25]). Let be a point at which singular and nonsingular sub arcs of an 
optimal control u are joined, and let q be the order of the singular arc. Suppose the strengthened GLC 
condition is satisfied at tc, and assume that the control is piecewise analytic in the neighborhood of tc- Let 
be the lowest order derivative ofu which is discontinuous at tc- Then q + r is an odd integer. 

In this problem, q is one because = H.y has the explicit formula with respect to 7 . The control 7 in 
the nonsingular arc has either 7 max or 7 min^ and hence 7 is discontinuous at the junction point, which means 
r = 0. Therefore, q + r = 1 in this problem. 

The singular arcs of B737-500 aircraft with various wind conditions are shown in Eig. 3. The right 
sides of singular arcs have negative values of rs(Vr, h) for both the fuel cost and the NOx cost cases. 
The singular controls of B737-500 with various wind conditions are shown in Eig. 4. As may be seen, 
discontinuity assumption at the junction point is true for any point on the singular arc. 


D. Boundary Arc 

In this subsection, we analyze the boundary arc, on which at least one of S{Vt, h) is zero. To do this, we 
assume that there exists a nonzero time interval [fi, t 2 ] in which the optimal trajectory is on the active state 
inequality. We also assume that the mixed state inequality C{Vt, h, 7 ) is not active on the boundary arc. 
Therefore, p{t) = 0,t £ [fi, f 2 ]- This assumption is reasonable because the flight path angle during the idle 
descent flight with maximum/minimum constant MACH/CAS is within the bound of the flight path angle 
limit. 
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Figure 3. B737-500 Singular arcs in various wind conditions 




Figure 4. B737-500 singular controls in various wind conditions 


From the necessary conditions for optimality and /r = 0, the following equations must be satisfied for 
the boundary arc interval: 


L^{VT,h,'y,Xv,Xh) = = 0 inEq. (39) 

Ho{VT,h, Xv, Xh) = 0 inEq. (42) 

SaiVr, h) = 0 
5a (Ft, f, 7) = 0 
: dH dSa 

^ ~ dVr 

• _ dH dSa 

^ dh dh 

7a(f) > 0, te[ti,t2] 

The active boundary arc S'a = 0 gives the relation between Vr and h. Since SaiVr, h) is the first order 
state inequality in this problem, which means that 7 related terms appear in the first time derivative Sa, 7 on 
the boundary arc can be obtained from the 5a = 0, and it has an explicit formulation as 


7 = - 


dSg 

dVr 


[g + VrWh,^) 


aSq 1"^ dSqP 
dh ^ dVr m 


(47) 
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The Lagrange multiplier of the active state inequality r/a is calculated by substituting the time derivative 
of \v obtained by differentiating Ay (Vr, h) into the adjoint equations of motion in Eq (46) 


OH , dSa 


-1 


la — ~(Ay + „ ) 


dVx'dVr 
m-f dSa 


-1 


> 0 


(48) 


where rs(Vr, h) is in Eq. (44). However, rs(Vr, h) can have a nonzero value on the boundary arc since 
rs(VT, h) = 0 only on the interior singular arc. 

In Eq. (48), E* (Vr, h) must be nonnegative to satisfy the necessary condition in Eq. (46) because 
is always positive during descent. By this analysis, the possible boundary arcs that satisfy the 
necessary condition are shown in Eig. 5. The bold lines represent the possible boundary arcs that can be a 
part of optimal trajectory when the singular arc crosses the state inequality bound. Though the CAS bound 
case in Eq. (11) is considered only in Eig. 5, the MACH bound case can be analyzed in the same way. 


I i 



01 

■a 


< 


ir.v = 0 

I 


dVr 


s,=o 

dS, 


<0 


S, =0 

dS, 


dv, 


>0 


•^1 “ Kii/i ^CAX 

5*2 — q — Lav 


\r,<o 

rv>o 


VcAs(Kt) 


■» 


Figure 5. Constrained arc satisfied necessary condition 


1. Continuity of adjoint variables 

The adjoint variables Ay and can be discontinuous for the boundary arc interval [ti,t 2 ] from the jump 
condition in Eq. (35) ~ (37). The adjoint variables can have a discontinuity at any time along the boundary 
arc interval[23]. Erom Eq. (33), Xx is constant along the whole optimal trajectory. Eurthermore, Ay and Xh 
are functions of Vr and h in Eq. (43) because and Hq are zero for the boundary arc interval, and hence 
these are also continuous along the boundary arc because the state variables are continuous. Therefore, the 
only possible discontinuous points are the junction points. The junction points have two types: the points 
at which the nonsingular arc and boundary arc are joined and the points at which the singular arc and the 
boundary arc are joined. To analyze the discontinuity of the adjoint variable at the junction points, we use 
the following two junction theorems in [26]. 

Theorem 2 (Theorem 5.1 in [26]). Let ti be the time at which an interior nonsingular arc and a boundary 
arc of an optimal control u are joined. Let u^'^\ r > 0, be the lowest order derivative of u which is 
discontinuous at ti and let p be the order of state inequality and q be the order of singular arc. Letp < 2q+r. 
Ifvfi) > 0, then p + r is an even integer. 
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Theorem 3 (Theorem 5.4 in [26]). Let ti be a point where an interior singular arc and a boundary arc of 
an optimal control u are joined. Let q be the order of the singular arc and assume that the strengthened 
GLC-condition holds. Let r > 0, be the lowest order derivative of u which is discontinuous at ti and 
let p < 2q + r. Then = 0, and q + r is an odd integer. 

In all the state inequalities in Eq. (11) and (12), p = 1 and the order of singular are q = 1. Therefore 
p < 2q + r holds always in this paper. On the nonsingular arc, 7 is either 7min or 7 max> and 7 stays in 
the interior of 0 set on the boundary arc by the assumption. Hence, 7 at the junction points, at which the 
nonsingular and the boundary arc are joined, is discontinuous. Therefore, r = 0 in this case, andp + r is odd 
number. According to Theorem 2, if p + r is odd number, then v{ti) = 0, which means adjoint variables are 
continuous at this type of junction points. According to Theorem 3, 1 / = 0 at the second type of junction 
points. Furthermore, r should be even number because q = 1 in this case. The singular arc in Fig. 3 shows 
that r = 0, and hence this condition is satisfied. By this analysis, we can conclude that adjoint variables are 
continuous along the optimal trajectory. 

IV. Optimal Trajectory Generation 

In this section, we present the optimal trajectory generation algorithm based on the analysis results in the 
previous section. In this algorithm, we generate the optimal trajectory by combining nonsingular, singular, 
and boundary arcs. The trajectory is generated by forward and backward integration without solving the 
optimal control problem. Therefore, it can be implemented in the FMS without additional computational 
capability. 

A. Algorithm 

The proposed optimal trajectory generation algorithm is as follows: 

STEP 1. Find the singular arc {Vt, h) using Eq. (44) (For this step, we solved the equation h) = 

0 numerically). 

STEP 2. Find the initial segment of the optimal trajectory by forward integrating the equations of motion 
in Eq. (1) ~ (3) from the initial point with either ymax or 7mm in LI{Vt, t). The control input is determined 
so that trajectory moves toward the singular arc. The integration stops when the trajectory reaches either the 
possible boundary arc or the singular arc. 

STEP 3. Find the last segment of the optimal trajectory by backward integrating the equations of motion 
from terminal point with either 'ymax or 7 mm in t). The control input is determined so that trajectory 

moves from singular arc to the terminal point. The backward integration stops when the trajectory reaches 
either the possible boundary arc or the singular arc. 

STEP 4. Determine the trajectory that is either on the possible boundary arc or the singular arc from 
section III.B and C. To do this, the trajectory is backward integrated from the junction point obtained from 
STEP 3 until it reaches the junction point from STEP 2. 

In STEP 4, if there exists a junction point between the boundary arc and singular arc, the trajectory 
changed from the singular arc to the boundary arc or vice versa depending on the possibility. 

The expected optimal trajectory structure is shown in Fig. 6 . The gray dash lines represent the possible 
boundary arcs by the optimality conditions in section III. The structure of the optimal trajectory depends on 
the initial and the terminal conditions as shown in Fig. 6 . Depending on the boundary conditions, boundary 
arc intervals may or may not exist. In case 1 of Fig. 6 , both junction points in STEP 2 and STEP 3 are on 
the boundary arcs, and hence the trajectory generated in STEP 4 consists of three segments: boundary arc, 
singular arc, boundary arc. On the other hand, the trajectory of STEP 4 is on the singular arc only in case 
2 of Fig. 6 . The number of segments and the type of segments differ from the boundary conditions. The 
structure of the optimal solution also depends on the wind profile because the singular arc is varied by the 
wind condition. 
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Figure 6. Switching structure of the optimal trajectory 


B. Optimality Check 

In this trajectory generation algorithm, we assume that the initial and the last parts of the trajectory, which 
moves toward the singular arc in STEP 2 and 3, are the parts of the optimal solution. To check if these parts 
are satisfied the necessary conditions, the sign of should be known. In this paper, is checked by the 
similar way in [14]. From the necessary condition, the adjoint variables Xy and are continuous at the 
junction point. Therefore, the values of Ay and at junction point are known from Eq. (28) and (42). From 
Eq. (34), Ay can be expressed as a function of Vr, h, and Xh- Then Eq. (25) is expressed as follows: 


‘-'■.(s'-'-S 


dKg,, ^ XhVrl + {Kdes - K,r{cVT + Wh)) dWhA 

dh + VrW/,^^) 7 ydh ^ dh j ' 


(49) 

From the continuity of Xh, the values of Xh at the junction points in STEP 2 and STEP 3 are known 
values. The value at the junction point in STEP 2 is the final value of the initial segment, and the value 
of STEP 3 is the initial value of the final segment. By integrating Eq. (49) with the known values at the 
junction points to the backward direction along the initial segment or to the forward direction along the 
final segmenf, Xh{t) can be calculated. From the relation between Ay and Xh, Ay also can be calculated. 
Therefore, the time history of H.y is obtained from Eq. (39), and the optimality condition can be checked by 
Eq. (38). The trajectory generated by the proposed algorithm satisfies this necessary condition, as will be 
shown by numerical examples. 


V. Numerical Example 

To verify that the trajectory generated by the proposed algorithm is the optimal, we compare it to the 
trajectory generated by the numerical method. For the numerical solution, we use the pseudospectral method 
since this method provided good solutions in the similar problems in [5]. GPOPS[27, 28] and SNOPT[29] 
are used as solvers for pseudospectral method and Nonlinear Programming (NEP) problem. 

We also analyze the wind effect on the optimal trajectory. The effects of the magnitude of the wind, 
wind shear, and the wind direction are evaluated. The numerical evaluations of the proposed algorithm are 
conducted with two aircraft types: B737-500 (B735) and B767-400 (B764). We use BADA[19] performance 
model for this numerical evaluation. However, since K^es and D are functions of Vt and h, which are quite 
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general, various performance models can be used to the formulation derived in section III. The boundary 
condition of this problem is given in Table 1. The maximum range dmax is set to -150 NM, which is large 
enough to satisfy the assumption in Lemma 1. The numerical values of the path constraints in Eq. (11) ~ 
(14) are given in Table 2. 


Table 1. Boundary condition for numerical example 


Initial Condition 

Final Condition 

VcAS (Kt) 

h (ft) (NM) 

ycAs (Kt) 

h(ft) 

(NM) 

265 

35,000 1 dmax 

250 

13,000 

-40 


Table 2. Path constraints of B735 and B764 


B737-500 

min max 

B767-400 

min max 

MACH 

0.45 

0.82 

0.45 

0.84 

CAS (Kt) 

220 

340 

230 

360 

ROD (m/s) 

-25.0 

-2.54 

-25.0 

-2.54 

7 (deg) 

-6 

0 

-6 

0 


A. Fuel Optimal Trajectory 

A comparison of the result with the proposed algorithm and with the numerical method is shown in Fig. 7. 
Zero wind is assumed in this comparison, however, we compared the results for various wind conditions, 
and the results for all the cases are similar. The two results in Fig. 7 are almost the same except for a 
very short interval near junction points calculated in STEP 1 and STEP 2 in the proposed algorithm. In the 
numerical method case, ^ 7 ^ term is added to the Eagrangian cost term in Eq. (19) to prevent the chattering 
phenomenon that usually occurs when the singular arc exists in the optimal trajectory [30]. This comparison 
shows that the proposed algorithm and the pseudospectral method generate the same trajectory. We checked 
the sign of using Eq. (49), and the results are described in Fig. 8 . In both initial and the final segments, 
the control inputs are ^max in Fig- 7. Therefore, should be negative to satisfy the necessary condition. 
The results show that the trajectory from the proposed method satisfies the necessary condition. 




a) Altitude profile b) CAS profile 

Figure 7. Comparison of numerical result and proposed algorithm result (fuel optimal case) 
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along track distance (NM) 

a) Initial segment 



along track distance (NM) 

b) Final segment 


Figure 8. H.y check for STEP 2 and STEP 3 in proposed algorithm 


1. Wind effect 

The fuel optimal trajectories of B735 and B764 for various wind speeds are shown in Fig. 9 and 10. The 
trajectory performances of the two aircraft are presented in Table 3 and 4. In these examples, wind was 
assumed to be a constant along track wind to evaluate wind speed effect. In both aircraft, the optimal TOD 
moves farther from runway threshold as wind speed increases. The descent speed during the singular arc 
decrease as the tail wind speed increases. 




Eigure 9. B735 fuel optimal trajectories with various wind conditions 


Table 3. B735 fuel optimal trajectory performance 


Wind(m/s) 

TOD(NM) 

TA(s) 

Fuel(kg) 

-30 

-96.441 

1154.489 

411.314 

-10 

-104.336 

1073.759 

342.119 

0 

-108.369 

1038.248 

311.588 

10 

-112.459 

1004.226 

283.161 

30 

-120.711 

942.495 

232.326 
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a) Altitude profile b) CAS profile 

Figure 10. B764 fuel optimal trajectories with various wind conditions 


Table 4. B764 fuel optimal trajectory performance 


Wind(m/s) 

TOD(NM) 

TA(s) 

Fuel(kg) 

-30 

-103.700 

1139.405 

729.267 

-10 

-111.959 

1057.429 

590.346 

0 

-116.167 

1020.959 

528.991 

10 

-120.419 

986.985 

472.219 

30 

-129.029 

925.355 

370.504 


2. Wind shear effect 

The fuel optimal CAS profiles of B735 for various wind profiles are shown in Fig. 11. The left column of 
Fig. 11 is the tail wind case and right column is the head wind case. To evaluate the effects of wind shear, 
the wind speeds at the cruise altitude are the same for all wind profiles. The influence of the wind shear on 
the optimal CAS profiles is relatively small. The results show that TOD points are almost the same for the 
various wind shear cases if the wind speeds at the cruise altitude are the same. 

3. Cross wind effect 

We analyzed the wind direction effect on the performance of the fuel optimal trajectories. For this analysis, 
we used two constant wind speed, 10 m/s and 30 m/s. The performance difference of B735 and B764 due 
to the wind directions are shown in Fig. 12. Time of Arrival (TA) and fuel burn increase as the along track 
wind speed decreases in both B735 and B764 cases. 

In order to compare the difference in optimal CAS profile, we compared the two wind profile cases. 
The first wind profile has the only along track component, and the second wind profile has the same along 
track wind speed as the first case and has a nonzero cross wind speed. The fuel optimal CAS profile for 
different cross wind profiles are shown in Fig. 13. In the case where the along track wind speed is zero, the 
differences in the optimal CAS profile between two cross wind speed cases are small, yet the differences 
are large in the case where the along track wind speeds are both 30 m/s, but the cross wind components are 
different (one is zero cross wind, and another has 60 m/s in magnitude, hence the angle ipui = 60 deg). In 
the second comparison case, which Fig. 13(b), the performance differences for the B735 are 15 sec in TA 
and 15 kg in fuel burn, while for the B764, they are 15 sec in TA and 35 kg in fuel burn. This result shows 
that we must take the cross wind component into account if we are to determine the true optimal trajectory. 
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Figure 11. Wind shear effects: left column - tail wind case, right column - head wind case 


Fuel burn (Kg) : B735 -d- lOin/s 

90 O 30ni/s 



Fuel burn (Kg) : B764 □ xom/s 


90 O 30m/s 



a) B735 Fuel 


b) B764 Fuel 


Figure 12. Cross wind effect on fuel optimal trajectories 
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a) B764, Wh '■ Om/s 



b) B764, Wh : 30m/s 


Figure 13. Fuel optimal CAS profiles with same along track winds 


B. Emissions Optimal Trajectory 

A comparison of the trajectory from the proposed algorithm and the optimal trajectory from the numerical 
method for minimizing NOx emissions is presented in Fig. 14. In this example, wind is assumed to be zero. 
The two trajectories are almost the same except small deviation at the junction point between nonsingular 
arc and singular arc. As in the fuel optimal case, the term is added with k = 0.01 to protect chattering 
for the numerical solution. It makes the solution smooth at the junction points compared to the solution 
from the proposed method. Since the NOx costs of the two trajectories are the almost the same, and the 
trajectories are almost the same, this smoothing effect at the junction points is negligible. Hence, it can be 
concluded that the proposed algorithm and the psuedospectral method generate the same trajectory. 

The minimum NOx trajectories for the B735 and B764 at various wind speeds are shown in Fig. 15 
and 16, respectively. The corresponding performance measures are presented in Table 5 and 6. To evaluate 
wind speed effect, we assume that there is a constant along track wind. The minimum NOx CAS profiles 
are different from the fuel optimal speed profiles. The CAS speeds are almosf consfanf along fhe singular 
arcs. The TOD poinfs are a liffle farther fhan fhaf for fhe fuel opfimal cases, and fhe arrival limes are slighlly 
increased compared lo fhe fuel opfimal cases. 

We analyzed fhe wind direcfion effecl on fhe minimum NOx frajeclory. The lesl condifions are fhe 
same as fhe minimum fuel frajeclory case. The performances of fhe opfimal solutions wifh fhe various wind 
direclions are shown in Fig. 17. As maybe seen, fhe NOx emission increases as fhe angle ijjw in Fig. 2 
increases. 

The cross wind effecl on fhe minimum NOx frajeclory is shown in Fig. 18. The comparison condifions 
are fhe same as fhe minimum fuel case in Fig. 13(b). The differences in fhe opfimal CAS profile of fhe 
Iwo wind profiles are large even fhough fhe along Irack wind speed of bofh cases are fhe same as 30m/s. 
Furlhermore, as shown in Fig. 18(a), fhe slruclure of fhe optimal speed profile can be differenl depending 
on fhe magnilude of fhe cross wind. The nonzero cross wind case has fhe boundary arc while fhe zero cross 
wind case does nof. This resulf furfher indicafes fhe imporfance of considering fhe cross wind lerm when 
deriving fhe opfimal frajeclory, regardless of objective. 
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Altitude (ft) Height (ft) 


Table 5. B735 minimum NOx trajectory performance 


Wind(m/s) 

TOD(NM) 

TA(s) 

iVOx(g) 

-30 

-97.008 

1195.025 

3698.059 

-10 

-104.858 

1107.863 

2932.396 

0 

-108.855 

1068.384 

2594.677 

10 

-112.890 

1031.335 

2282.597 

30 

-121.035 

960.993 

1725.159 


Table 6. B764 minimum NOx trajectory performance 


Wind(m/s) 

TOD(NM) 

TA(s) 

iVOx(g) 

-30 

-104.017 

1169.034 

6700.151 

-10 

-112.262 

1083.306 

5123.663 

0 

-116.457 

1044.782 

4428.391 

10 

-120.692 

1008.701 

3785.831 

30 

-129.261 

943.040 

2636.919 


xlO* 




a) Altitude profile b) CAS profile 

Figure 14. Comparison of numerical method and proposed algorithm results (NOx optimal case) 
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Figure 15. B735 minimum NOx trajectories with various wind conditions 
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Figure 16. B764 minimum NOx trajectories with various wind conditions 
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Figure 17. Cross wind effect on NO:c optimal trajectories 
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VI. Conclusion 


We analyzed the optimal vertieal trajeetories that minimize the environmental impaets of a CDA in the 
presenee of both along traek and eross wind eomponents. In this paper, both eross wind and along traek 
wind were assumed to be funetions of altitude. Flight idle thrust was assumed during the entire deseent 
phase. The flight range was speeified from a point during the latter stages of the eruise segment to the meter 
fix. We eonsidered CAS and Maeh eonstraints(whieh are the pure state path eonstraints), flight path angle 
eonstraints, and a maximum deseent rate limit(whieh is a mixed input and state path eonstraint). 

We formulated this problem as a single phase initial eondition free optimal eontrol problem. The for¬ 
mulation of the eost funetional that we used is quite general sinee it is a funetion of airspeed and altitude, 
and this formulation was used for the fuel eost and the gaseous emissions eost minimization problems. We 
found an explieit formulation of the singular by using the neeessary eondition for the optimality of the op¬ 
timal eontrol problem ineluding path eonstraints. We also found the existing eondition of the boundary are 
on whieh CAS/Maeh is eonstant. From these analyses, we proposed an algorithm to generate the optimal 
trajeetory that minimizes the eost funetional by forward and baekward integration. The optimal trajeetories 
generated by the proposed algorithm were evaluated by eomparing to the optimal trajeetories obtained from 
the numerieal method. Sinee the proposed algorithm does not require additional eomputational power rela¬ 
tive to the eurrent trajeetory generation method using the VNAV funetion in the FMS, this algorithm eould 
easily be implemented in the FMS for the online trajeetory optimization. 
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